Distinct H3K27me3 and H3K27ac Modifications in Neural Tube Defects Induced by Benzo[a]pyrene

The pathological mechanisms of neural tube defects (NTDs) are not yet fully understood. Although the dysregulation of histone modification in NTDs is recognized, it remains to be fully elucidated on a genome-wide level. We profiled genome-wide H3K27me3 and H3K27ac occupancy by CUT&Tag in neural tissues from ICR mouse embryos with benzo[a]pyrene (BaP)-induced NTDs (250 mg kg−1) at E9.5. Furthermore, we performed RNA sequencing (RNA-seq) to investigate the regulation of histone modifications on gene expressions. Gene ontology and KEGG analysis were conducted to predict pathways involved in the development of NTDs. Our analysis of histone 3 lysine 27 modification in BaP-NTD neural tissues compared to BaP-nonNTD revealed 6045 differentially trimethylated regions and 3104 acetylated regions throughout the genome, respectively. The functional analysis identified a number of pathways uniquely enriched for BaP-NTD embryos, including known neurodevelopment related pathways such as anterior/posterior pattern specification, ephrin receptor signaling pathway, neuron migration and neuron differentiation. RNA-seq identified 423 differentially expressed genes (DEGs) between BaP-NTD and BaP-nonNTD group. The combination analysis of CUT&Tag and RNA-seq found that 55 DEGs were modified by H3K27me3 and 25 by H3K27ac in BaP-NTD, respectively. In the transcriptional regulatory network, transcriptional factors including Srsf1, Ume6, Zbtb7b, and Cad were predicated to be involved in gene expression regulation. In conclusion, our results provide an overview of histone modifications during neural tube closure and demonstrate a key role of genome-wide alterations in H3K27me3 and H3K27ac in NTDs corresponding with changes in transcription profiles.


Introduction
Neural tube defects (NTDs) are among the most common and serious birth defects, caused by the failure of the neural tube closure between 21 and 28 days after conception [1]. According to the World Health Organization, approximately 300,000 babies are born each year with NTDs [2], resulting in approximately 88,000 deaths and 8.6 million disabilityadjusted life years [3,4]. In low-income countries, NTDs may account for 29% of neonatal

Animal Experiments
The BaP-induced NTD mouse model was performed as previously described [10]. Briefly, ICR mice of 8-9 weeks old weighing 28 ± 2 g were used in the experiment. Female mice were mated with males overnight and vaginal plugs were examined the following morning. Noon on the day of finding a vaginal plug was considered 0.5 days of embryonic development (E0.5). Pregnant mice were randomly divided into 2 groups. In BaP-treatment groups, mice were given BaP (CAS-No. 192-97-2, purity 98%, Sigma-Aldrich, St. Louis, MO, USA) intraperitoneally, dissolved in corn oil (CAS-No. 8001-30-7, aladdin, Shanghai, China) (25 mg mL −1 ), from E6.5 for three consecutive days at a dose of 250 mg kg −1 . Mice in the control group were treated with corn oil (10 mL kg −1 ). On E9.5, pregnant mice were sacrificed by cervical dislocation and the embryos were removed by caesarean section. Embryos were carefully inspected for visible external malformations under a dissecting microscope. NTD-affected embryos were classified as showing distinct evidence of failed closure of the cephalic neural tube. Finally, 16.4% of the embryos were observed with open neural tube (Table S1). The somite numbers were carefully counted. E9.5 embryos with 22-27 somites were used for further analysis [19].

Isolation of Neural Tissues
In BaP-treatment group, E9.5 embryos affected with cephalic NTDs were collected as BaP-NTD group, and embryos from the same dam without obvious malformations were collected as BaP-nonNTD group. Embryos with no malformations from the control group were collected as control group. Specially, neural tissues of five E9.5 embryos collected from five independent dams were pooled as one sample for RNA-seq experiment, and neural tissues from fourteen embryos collected from seven independent dams were pooled as one sample for CUT&Tag experiment. Embryonic neural tissues, from the most rostral aspect of the forebrain to the caudal aspect of the hindbrain, were excised according to the previous study [20] and further trimmed to eliminate any non-neural tissues as precisely as possible, the precise scope was shown in Figure 1. Three biological replicates for each group were conducted for RNA-seq and CUT&Tag, respectively.
showing distinct evidence of failed closure of the cephalic neural tube. Finally, 16.4% of the embryos were observed with open neural tube (Table S1). The somite numbers were carefully counted. E9.5 embryos with 22-27 somites were used for further analysis [19].

Isolation of Neural Tissues
In BaP-treatment group, E9.5 embryos affected with cephalic NTDs were collected as BaP-NTD group, and embryos from the same dam without obvious malformations were collected as BaP-nonNTD group. Embryos with no malformations from the control group were collected as control group. Specially, neural tissues of five E9.5 embryos collected from five independent dams were pooled as one sample for RNA-seq experiment, and neural tissues from fourteen embryos collected from seven independent dams were pooled as one sample for CUT&Tag experiment. Embryonic neural tissues, from the most rostral aspect of the forebrain to the caudal aspect of the hindbrain, were excised according to the previous study [20] and further trimmed to eliminate any non-neural tissues as precisely as possible, the precise scope was shown in Figure 1. Three biological replicates for each group were conducted for RNA-seq and CUT&Tag, respectively.

RNA-seq
Total RNA was extracted from neural tissues using the TRIZOL method and the concentration was determined using NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA). Library preparation was performed using the Illumina's NEBNext ® UltraTM RNA Library Prep Kit (Illumina, San Diego, CA, USA) and high-throughput RNA-seq was performed using the Illumina Novaseq 6000 platform (Illumina, San Diego, CA, USA) with 100 bp paired-end sequencing.

Analysis of Differential Gene Expression and Functional Annotation
RNA-seq raw reads were removed adapters and trimmed low quality reads with Trimmomatic v.0.36. These processed reads were then mapped to the mouse genome (GRCm38/mm 10) using Hisat2 software. The readings were counted using HTseq-count and normalized to fragments per kilobase (FPKM) for further visualization [21]. The

RNA-seq
Total RNA was extracted from neural tissues using the TRIZOL method and the concentration was determined using NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA). Library preparation was performed using the Illumina's NEBNext ® UltraTM RNA Library Prep Kit (Illumina, San Diego, CA, USA) and high-throughput RNA-seq was performed using the Illumina Novaseq 6000 platform (Illumina, San Diego, CA, USA) with 100 bp paired-end sequencing.

Analysis of Differential Gene Expression and Functional Annotation
RNA-seq raw reads were removed adapters and trimmed low quality reads with Trimmomatic v.0.36. These processed reads were then mapped to the mouse genome (GRCm38/mm 10) using Hisat2 software. The readings were counted using HTseq-count and normalized to fragments per kilobase (FPKM) for further visualization [21]. The differentially expressed genes (DEGs) were identified using the DESeq2 package and Brain Sci. 2023, 13, 334 4 of 16 defined as a fold change (FC) >2 and a p-value < 0.05. The DEGs were then subjected to Gene Ontology (GO) and pathway enrichment analysis to identify the significantly enriched functional classification, signaling and metabolic pathways, which were performed based on Gene Ontology Database (http://www.geneontology.org/, accessed on 1 January 2019) [22] and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (http://www.genome.jp/kegg/, accessed on 1 June 2022) [23] with R package clusterProfiler. GO terms were identified to be significantly enriched when p < 0.05 after Benjamini multiple test correction. Pathways were identified to be significantly enriched when false discovery rate (FDR) < 0.05.

CUT&Tag Data Processing and Analysis
We used FastQC to check the quality of the raw reads. Low quality reads as more than 10% of bases with quality scores less than 20 were filtered. The remaining reads that passed all the filtering steps were counted as clean reads and aligned to the mouse genome (GRCm38/mm 10 build) with the BWA mem v.0.7.12. Only non-redundant and uniquely mapped reads were retained to correct for sequencing bias. Peaks and enriched regions were called with MACS2 v.2.1.0 [24]. Noisy peaks with very weak signals were removed in further analysis. PeakAnnotator was used to identify the nearest TSS of every peak and the distance distribution between peaks and TSS was shown. Moreover, the distribution of peak summits on different function regions, such as 5utr, CDS, 3utr, was performed. Different peak analysis was based on the software PePr v.1.1.18 and defined as the FC between two groups being more than 2, p < 1 × 10 −5 [24]. Peak-related genes were confirmed by PeakAnnotator. Data were visualized in the UCSC genome browser. Functional analyses were performed using the edgeR package. Hypergeometric optimization of motif enrichment (HOMER) was used to discover motifs of specific regions [25].

Quantitative Real-Time Quantitative Polymerase Chain Reaction (qRT-PCR)
RNA was isolated from neural tissues of E9.5 embryos using Trizol (Invitrogen, Carlsbad, CA, USA); genomic DNA was removed by DNase I digestion (DNA-free, Invitrogen, Carlsbad, CA, USA) and then reverse-transcribed using random hexamers (Superscript VILO cDNA synthesis kit, Invitrogen, Carlsbad, CA, USA). The abundance of mRNA of Myl2, Csrp3 and Grm2 were analysed using real-time PCR (iTaqTM Universal SYBR Green Supermix, BioRad, Hercules, CA, USA) on a 7500 Fast Real Time PCR system (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, USA), with each sample analysed in triplicate. Primers are listed in Table S2. Relative quantification of each gene expression level was normalized according to the Gapdh gene expression.

Statistical Analyses
Statistical analyses were carried out using SPSS version 20.0 (IBM, Chicago, IL, USA). Non-paired Student's t-test was used when two groups were being compared. Differences were considered statistically significant at p < 0.05. Data are represented as the mean ± standard error of mean (SEM) of at least in triplicate.

CUT&Tag Analysis
Neural tissues from three experimental groups: NTDs in BaP exposed group (BaP-NTD), non-NTDs in BaP exposed group (BaP-nonNTD), and normal embryos in control group (control) were extracted to profile epigenetic changes ( Figure 1a). The quality of the clean reads and mapping rates of CUT&Tag data are shown in Tables S3 and S4. Figure 1b-d summarized differences in H3K27me3 and H3K27ac abundance between the three different experimental groups. BaP exposure significantly decreased the total number of H3K27me3 peaks, while had little impact on the number of H3K27ac peaks. The H3K27me3 and H3K27ac peak widths in BaP-NTD group were wider than those in control and BaP-nonNTD group (Table S5). The distribution of genomic regions modified by H3K27me3 and H3K27ac was classified into eight regions (promoters, upstream2k, 5utr, exons, introns, 3utr, downstream2k and cgi) according to the UCSC Genome Browser. The H3K27me3 and H3K27ac marks were mainly mapped to intron and promoter regions ( Figure 1d). BaP exposure resulted in decreased H3K27me3 peaks in upstream2k, intron and cig, while it had little to no impact on H3K27ac peaks across the whole gene ( Figures S1 and S2). No significant differences on H3K27me3 or H3K27ac peak distributions were observed between BaP-NTD and BaP-nonNTD embryos.
As shown in Figure 2a,b, there were unique H3K27me3 binding genes present in each group. Compared to controls, 397 mutual H3K27me3 and 226 H3K27ac binding genes were identified in both BaP-NTD and BaP-nonNTD group, and 895 H3K27me3 and 431 H3K27ac binding genes were unique for BaP-NTD group. The GO enrichment analysis showed that these 397 mutual H3K27me3 peak related genes were mainly involved in synapse organization, axonogensis, dendrite development and postsynaptic membrane ( Figure 2c). Similar enrichment terms were observed for the 226 mutual H3K27ac binding genes ( Figure 2d, Tables S6 and S7). Function annotation of the 895 H3K27me3 binding genes uniquely associated with BaP-NTD group revealed a significant enrichment of terms associated with neurodevelopment including growth cone, site of polarized growth and distal axon (Figure 2e, Table S8). However, no significant enrichment pathways were observed for these 431 H3K27ac unique binding genes ( Figure 2f, Table S9). The top 10 unique peak binding genes identified in BaP-NTD are listed in Figure 2g,h.
A total of 6879 and 9914 differential H3K27me3 peaks, and 1574 and 3753 differential H3K27ac peaks were identified in BaP-NTD/control and BaP-nonNTD/control comparison, respectively (Tables S10-S13). The overview of the enriched pathways based on these differential peak-associated genes is presented in Figure 3a. Compared to BaP-nonNTD/control comparison, more pathways associated with neural development were enriched for BaP-NTD/control comparison. When compared the histone modification between BaP-NTD and BaP-nonNTD, 6045 differential H3K27me3 peaks and 3104 differential H3K27ac peaks were identified (Tables S14 and S15). As shown in Figure 3c,d, nerve system development was enriched as one of the top terms for both the differential H3K27me3 and H3K27ac binding genes for BaP-NTD/BaP-nonNTD comparison. Other classical neurodevelopment related pathways enriched included neurogenesis, neuron migration, neuron differentiation, neural crest cell migration, axon guidance, axonogenesis, brain development, brain morphogenesis, forebrain development, cerebellum development, cerebral cortex development, anterior/posterior pattern specification, ephrin receptor signaling pathway (Tables S16 and S17). Representative differential H3K27me3 and H3K27ac peaks between BaP-NTD and BaP-nonNTD are shown in Figure 3b A total of 6879 and 9914 differential H3K27me3 peaks, and 1574 and 3753 differential H3K27ac peaks were identified in BaP-NTD/control and BaP-nonNTD/control comparison, respectively (Tables S10-S13). The overview of the enriched pathways based on these differential peak-associated genes is presented in Figure 3a. Compared to Top six GO biological processes, cellular components, and molecular functions enriched from common H3K27me3 (c) and H3K27ac (d) binding genes between BaP-nonNTD and BaP-NTD (ordered by significance level). Top six GO biological processes, cellular components, and molecular functions enriched from unique H3K27me3 (e) and H3K27ac (f) binding genes for BaP-NTD (ordered by significance level). The list of the top 10 unique BaP-NTD peak binding genes for H3K27me3 (g) and H3K27ac (h). axonogenesis, brain development, brain morphogenesis, forebrain development, cerebellum development, cerebral cortex development, anterior/posterior pattern specification, ephrin receptor signaling pathway (Tables S16 and S17). Representative differential H3K27me3 and H3K27ac peaks between BaP-NTD and BaP-nonNTD are shown in Figure 3b. Integrative genomics viewer snapshot of representative differential modified H3K27me3 and H3K27ac signals in comparison of BaP-NTD with BaP-nonNTD. Biological process GO enrichment terms for the differential H3K27me3 (c) and H3K27ac (d) peak binding genes identified in BaP-NTD-vs-BaP-nonNTD comparisons, ordered by adjusted p-value. Neurodevelopment associated terms were marked with red circles.

Transcriptome Analysis
To further investigate how H3K27 modification status links with individual gene expression, we created a total RNA-seq data. An overview of the RNA-seq data is presented in Table S18. The relationship between histone modifications and gene expressions across all the three groups is presented in Figure 4a. As expected, genes with an H3K27me3 mark in the promoter region had a lower expression, those with an H3K27ac mark had a higher expression, and those with both marks had a distribution of expression levels somewhere in between.
A total of 1501 DEGs were observed for BaP-NTD/control comparison, which was significantly higher than that for BaP-nonNTD/control comparison (673 DEGs) ( Figure   Figure 3. Analysis of differentially modified histone peaks in E9.5 mouse embryos with and without NTDs. (a) Overview of KEGG enrichment analysis of differentially modified histone peak binding genes identified in BaP-NTD-vs-control and BaP-nonNTD-vs-control comparisons. (b) Integrative genomics viewer snapshot of representative differential modified H3K27me3 and H3K27ac signals in comparison of BaP-NTD with BaP-nonNTD. Biological process GO enrichment terms for the differential H3K27me3 (c) and H3K27ac (d) peak binding genes identified in BaP-NTD-vs-BaP-nonNTD comparisons, ordered by adjusted p-value. Neurodevelopment associated terms were marked with red circles.

Transcriptome Analysis
To further investigate how H3K27 modification status links with individual gene expression, we created a total RNA-seq data. An overview of the RNA-seq data is presented in Table S18. The relationship between histone modifications and gene expressions across all the three groups is presented in Figure 4a. As expected, genes with an H3K27me3 mark in the promoter region had a lower expression, those with an H3K27ac mark had a higher expression, and those with both marks had a distribution of expression levels somewhere in between.
A total of 1501 DEGs were observed for BaP-NTD/control comparison, which was significantly higher than that for BaP-nonNTD/control comparison (673 DEGs) (Figure 4b and Tables S19 and S20). The expressions of five representative neural tube closure related genes among the three groups are shown in Figure 4c. Compared to the control group, the expressions of Frem2, Vangl1, Ptk7, were significantly increased, while the expressions of Pax7, and Fkbp8 were decreased in the BaP-NTD group. In contrast, the levels of these genes in the BaP-nonNTD group were somewhere in between, and 423 DEGs were found for BaP-NTD/BaP-nonNTD comparison (Figure 4b, Table S21). Over-all, biological function analysis revealed that these 423 DEGs were mainly classified as three classes: (1) genes involved in the cellular process/metabolic process, mainly for neurodevelopment (314 genes), (2) genes involved in multicellular organismal process (206 genes), and one of the top terms was nervous system development (3) genes involved in localization (153 genes) (Figure 4d,e).
group, the expressions of Frem2, Vangl1, Ptk7, were significantly increased, while the expressions of Pax7, and Fkbp8 were decreased in the BaP-NTD group. In contrast, the levels of these genes in the BaP-nonNTD group were somewhere in between, and 423 DEGs were found for BaP-NTD/BaP-nonNTD comparison (Figure 4b, Table S21). Overall, biological function analysis revealed that these 423 DEGs were mainly classified as three classes: (1) genes involved in the cellular process/metabolic process, mainly for neurodevelopment (314 genes), (2) genes involved in multicellular organismal process (206 genes), and one of the top terms was nervous system development (3) genes involved in localization (153 genes) (Figure 4d,e). Gene ontology analysis of common differentially expressed genes identified in BaP-NTD-vs-BaP-nonNTD comparisons. GO terms associated with cellular process/metabolic process (yellow), multicellular organismal process (blue), and localization (grey) were highlighted respectively. (e) Genes contributing to the three main GO term classes were divided into up-(red) and down-regulated (gray) fractions. Gene numbers in each fraction are indicated.

Combined Analysis of Histone Modifications and Transcription
If we further characterize the DEGs with differential H3K27 modifications, no apparent trends were observed between gene expressions with histone modifications (Figures 5a, S3 and S4). A total of 198 DEGs with H3K27me3 modification and 65 DEGs with H3K27ac modification unique for BaP-NTD/control comparison were observed ( Figure  5b, Tables S22 and S23). The expressions of three top genes, Myl2, Csrp3 and Grm2,

Combined Analysis of Histone Modifications and Transcription
If we further characterize the DEGs with differential H3K27 modifications, no apparent trends were observed between gene expressions with histone modifications (Figures 5a, S3 and S4). A total of 198 DEGs with H3K27me3 modification and 65 DEGs with H3K27ac modification unique for BaP-NTD/control comparison were observed (Figure 5b, Tables S22 and S23). The expressions of three top genes, Myl2, Csrp3 and Grm2, were validated by qRT-PCR ( Figure S5). When compared BaP-NTD to BaP-nonNTD, 55 DEGs with H3K27me3 mod-ification and 25 DEGs with H3K27ac modification were observed (Table S24). We then performed GO analysis to investigate the biological function and discovered that these 55 DEGs with peaks of H3K27me3 were highly enriched in growth and developmental growth, and these 25 DEGs with H3K27ac modification were enriched in growth and neuron differentiation (Figure 5c). Combining the 263 DEGs unique for BaP-NTD/control comparison with the 80 DEGs for BaP-NTD/BaP-nonNTD comparison, 19 common DEGs including Myl2, Dlg2 and Rgs6 were observed ( Table 1).

of 16
and discovered that these 55 DEGs with peaks of H3K27me3 were highly enriched in growth and developmental growth, and these 25 DEGs with H3K27ac modification were enriched in growth and neuron differentiation (Figure 5c). Combining the 263 DEGs unique for BaP-NTD/control comparison with the 80 DEGs for BaP-NTD/BaP-nonNTD comparison, 19 common DEGs including Myl2, Dlg2 and Rgs6 were observed (Table 1).

Motif Analysis
We further performed motif analysis to identify the potential transcription factors which may regulate the DEGs in BaP-NTDs by HOMER, the de novo motif analysis. As shown in Figure 5d, five motifs were significantly (threshold, p-value < 1 × 10 −12 ) enriched in upregulated genes associated with increased H3K27me3 peaks, and two motifs in upregulated genes with decreased H3K27me3 peaks for BaP-NTD/control comparison. Of these, Smad2 and SRSF1 were previously shown to be involved in the development of NTDs. No significant transcription factor motifs were identified in other groups (Tables S25-S27).

Discussion
In this study, we correlated histone modifications to gene expressions in mouse NTDs via integrative analysis of the RNA-seq and CUT&Tag data, with both internal controls and vehicle controls to provide new and comprehensive information on epigenetic regulation during neural tube closure. Our data indicated several genes with specific alterations of histone modifications in NTDs that might be responsible for the failure of neural tube closure, such as Myl2, Dlg2, and Rgs6. Function analysis revealed potential key pathways implicated in the pathological process of NTDs, such as anterior/posterior pattern specification, ephrin receptor signaling pathway, neuron migration and neuron differentiation. In addition, our analysis identified new transcription factors which may be involved in the epigenetic regulation of molecular events in BaP-induced NTDs, including Srsf1, Ume6, Zbtb7b, and Cad.
Previous findings suggested that PAHs exposure was associated with the occurrence of NTDs [11,26]. However, the underlying pathological mechanisms have not been fully elucidated. Neurulation is a complex and multistep process, involving the precise temporal and spatial regulation of gene expression [27]. A considerable number of biological events are involved in the process of neural tube closure, such as apoptosis, proliferation, differentiation and migration [27,28]. Thus, it is of great significance to explore the key molecular events and regulatory mechanisms that lead to the failure of neural tube closure. With the comparison of internal and vehicle control group, our results identified a number of pathways specially enriched in NTDs, such as anterior/posterior pattern specification, which was closely related with neurodevelopment [29]. Moreover, we found that the ephrin receptor signaling pathway might be implicated in BaP induced NTD. Ephrin-A5 is a recognized gene implicated in cranial NTD based on mice models [30]. Evidence was further presented that Ephrin/Eph-mediated signaling was tightly associated with neural tube anterior-posterior patterning [31]. However, this research was performed with vertebrate spinal cord. Whether Ephrin participated in anterior-posterior patterning during cranial neural tube closure required further studies.
In addition, the GO analysis in the present study revealed that compared to its internal controls, a number of pathways participated in neuron differentiation and migration were mainly enriched in embryos affected by NTDs, especially in the forebrain. Our previous work has shown that the majority of NTDs induced by BaP exposure occurred in the forebrain [10]. During embryonic development, the tightly regulated cellular events including proliferation, apoptosis, differentiation, and migration are crucial for neural tube development [32,33]. The involvement of neuron apoptosis in BaP-induced NTD was confirmed in our previous study [10]. Data presented here indicated that the disruption of neuron differentiation and migration might be also one of the main underlying mechanisms for neural tube closure failure caused by BaP. Recent studies have confirmed the effect of histone medication on neurulation through regulation of cell differentiation. One study from Patterson et al. found that in Mox2-Cre mutant embryos, proliferation and apoptosis was not affected in neural tissues, while the disruption of interneuron differentiation and specification was the major contributor to the observed NTDs, and they further suggested the involvement of histone acetylation in the regulation of neuronal differentiation [34]. Li et al. reported that the EZH2-mediated H3K27 methylation was responsible for the premature differentiation of neuron cells in zebrafish with deficits in forebrain caused by knockdown of Sox19b [35]. Another NTD zebrafish model induced by knockdown of Pcgf1 also indicated the role of histone modification in regulation of neuron differentiation [36]. However, relevant studies on the epigenetic regulation of neuron migration through histone modification in NTD were scarce. Further research is needed to fully clarify the underlying mechanisms of dysregulation of neuron differentiation and migration caused by histone modification in NTD.
Deciphering the molecular underpinnings of neural tube development would be helpful for providing new strategies for prevention. We identified a set of genes with changes of both transcriptions and histone modifications specifically in BaP-NTD group, some of which were reported to be implicated in neurological disorders. For example, Ctnna3 was implicated in Alzheimer's disease [37]. Inhibition of Dpp4 showed a neuroprotection effect in diabetic mouse brain [38], and the Cacna2d2 mutant mice were affected with cerebellar neuro-degeneration associated with ataxia and seizures [39]. However, the role of these genes in neurodevelopment during embryonic stage was not clear. Moreover, we also found several candidate genes that were associated with neuron development. For example, our data showed an increased expression of Dlg2 with decreased of H3K27me3 specifically in BaP-NTD group, which was a known gene vital for neuronal differentiation, maturation, and migration [40]. Myl2 was reported as one of the main down regulated genes in In Adnp −/− NTD mice, and Rgs6 was one of the upstream modulators of Notch signaling, and the mutant of which would result malformations including delayed neural tube maturation [41]. The present work indicated higher expression of Myl2 and Rgs6 with histone modification in BaP-NTD group compared to BaP-NTD and control group. Thus, our results implied that the up-regulation of these two genes might also cause NTD. Altogether, our work provided new evidence for the regulation role of H3K27 histone modifications in the etiology of NTDs and suggested potential new gene targets for NTDs pathogenesis.
Interestingly, by de novo motif analysis, we identified enrichment of sequence motifs unique for the differentially trimethylated regions in NTD embryos. Among them, Smad2 has been well demonstrated to be vital for neural tube closure by regulating neural epithelial organization, including tissue shape, size and three-dimensional morphogenesis through tight junction pathways [42]. In frog, Smad10 was required for the formation of nervous system, as a key component in the development of both anterior and posterior neural tissue [43]. Srsf1 has emerged as a key oncodriver in tumor progression such as gliomagenesis for its role in alternative splicing [44]. The target genes of Srsf1 included Stat3, Bcl-x, Mdm2 [45][46][47], which were all related to neural tube closure [48][49][50]. Moreover, Ume6 was used to be regarded as a key transcriptional regulator of early meiotic gene expression [51,52], and its role in autophagy has been documented recently [53]. Zbtb7b belongs to the large ZBTB family, known as critical factors regulating the lineage commitment and differentiation of hematopoietic-derived cells [54]. Cad is a multifunctional protein that participates in the initial three speed-limiting steps of pyrimidine nucleotide synthesis. The important role of Cad in neurological disorders including epileptic encephalopathy and Alzheimer's disease has also been emphasized [55]. However, the direct molecular function of Ume6, Zbtb7b and Cad in neural tube closure has not been well characterized.
Few studies have examined the effects of BaP exposure on neurodevelopment in offspring without obvious NTDs. Air pollution continues to be a worldwide environmental health problem, which has an enormous influence on fetal development beginning during gestation. PAHs are one of the main toxic, mutagenic, and carcinogenic air pollutants and can cross the placenta. Our study found that although some embryos showed no obvious neural tube malformations after BaP treatment, their neurodevelopment were impaired as shown by the differentially expressed and histone modified neurodevelopment-related genes, such as wnt7b, Fzd and Cdh9. Epidemiological studies have demonstrated that PAHs exposure were positively correlated with the increased risks for neurobehavioral development problems, including attention deficit hyperactivity disorder [56], decrease of IQ [57], and cognitive developmental delay [58]. Animal experiments implied that the potential underlying mechanisms included oxidative stress, inflammation, altered levels of dopamine and/or glutamate, and changes in synaptic plasticity/structure [59]. The synapse is fundamental for the sophisticated ensemble of the highly specialized neuron network. The non-malformed BaP-treated embryos in the present study showed significant disturbances in pathways associated with synapse, including synapse organization, synapse assembly, chemical synaptic transmission and regulation of synaptic plasticity. Our data supported the previous results, emphasizing the role of synapse in PAHs-induced embryo neurotoxicity, and provided new insights for further study.
It should be mentioned that there are some limitations in this study. The candidate genes and related pathways found in our study were not further investigated. Functional analyses of key genes in further studies were required to obtain more in-depth understanding for development of NTDs. Another drawback is that we used intraperitoneal injection of BaP to induce NTDs in mice, while humans are mostly exposed to PAHs through inhalation and diet in the real world. Previous works from ours and others have tried to induce NTDs with BaP through gavage, however, few NTDs were observed [60,61]. We acknowledged that the pathogenesis of NTDs in mice induced by intraperitoneal injection of BaP could be different from that of human NTDs. However, many of the identified pathways and genes in the present study were reported to be implicated in other NTD animal models and human studies. Our previous studies have identified some key pathogenic genes and pathway shared between this BaP-NTD mouse model and human cases [10,62].

Conclusions
With this study, we showed that genome-wide H3K27me3 and H3K27ac modification profiles were useful to discover novel molecular pathways and genes for the better understanding of NTD pathophysiology. We would recommend future studies to confirm the newly identified targets in the present study. Moreover, since neural tissue samples represent an average state across different cellular compartments, newer techniques looking at subpopulations of cells might help to unravel cell-specific deregulated pathway.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/brainsci13020334/s1, Figure S1: Genome-wide peak distribution of H3K27me3 in E9.5 embryos with and without neural tube defects; Figure S2: Genome-wide peak distribution of H3K27ac in E9.5 embryos with and without neural tube defects; Figure S3: The relationships between fold change of CUT&Tag signals for H3K27me3/H3K27ac and fold change of gene expressions in BaP-NTD-vs-control comparison; Figure S4: The relationships between fold change of CUT&Tag signals for H3K27me3/H3K27ac and fold change of gene expressions in BaP-nonNTD-vs-control comparison; Figure S5: Relative expression of mRNA of E9.5 embryos E9.5 embryos with and without neural tube defects by real-time PCR. Table S1: Embryotoxicity of BaP in ICR mice; Table S2: Sequences of primers for real-time PCR; Table S3: Descriptive summary of data generated by CUT&Tag for H3K27me3; Table S4: Descriptive summary of data generated by CUT&Tag for H3K27ac; Table S5: Peak widths of H3K27me3 and H3K27ac in E9.5 mouse embryos with and without neural tube defects; Table S6: Gene Ontology enrichment analysis of 397 mutual H3K27me3 peak binding genes for BaP-NTD and BaP-nonNTD; Table S7: Gene Ontology enrichment analysis of 226 mutual H3K27ac peak binding genes for BaP-NTD and BaP-nonNTD; Table S8: Gene Ontology enrichment analysis of 895 unique H3K27me3 peak binding genes for BaP-NTD; Table S9: Gene Ontology enrichment analysis of 431 unique H3K27ac peak binding genes for BaP-NTD; Table S10: Lists of H3K27me3 differential peaks between BaP-NTD and control group; Table S11: Lists of H3K27me3 differential peaks of BaP-nonNTD and control group; Table S12: Lists of H3K27ac differential peaks of BaP-NTD and control group; Table S13: Lists of H3K27ac differential peaks BaP-nonNTD and control group; Table S14: Lists of H3K27me3 differential peaks between BaP-NTD and BaP-nonNTD group; Table S15: Lists of H3K27ac differential peaks between BaP-NTD and BaP-nonNTD group; Table S16: Gene Ontology enrichment analysis of differentially H3K27me3 peak binding genes for BaP-NTD compared with BaP-nonNTD; Table S17: Gene Ontology enrichment analysis of differentially H3K27ac peak binding genes for BaP-NTD compared with BaP-nonNTD; Table S18: Descriptive summary of data generated by RNA-seq; Table S19: Lists of DEGs between BaP-NTD and controls; Table S20: Lists of DEGs between BaP-nonNTD and controls; Table S21: Lists of DEGs between BaP-NTD and BaP-nonNTD; Table S22: Lists of the overlaps of differential expressed genes and histone modified genes of BaP-NTD vs. control group; Table S23: Lists of the overlaps of differential expressed genes and histone modified genes of BaP-nonNTD vs. control group; Table S24: Lists of the overlaps of differential expressed genes and histone modified genes of BaP-NTD vs. BaP-nonNTD; Table S25: Results of motif enrichment by HOMER of BaP-NTD vs. control group; Table S26: Results of motif enrichment by HOMER of BaP-nonNTD vs. control group; Table S27: Results of motif enrichment by HOMER of BaP-NTD vs. BaP-nonNTD group.

Conflicts of Interest:
The authors declare no conflict of interest.